function istp=sum_tstep(istp)
% sumstep has been re-written as a function. 
% input = time step you want to watch
%
%
istp=2200;
xmin = -30;
xmax = 30;

zmin = -20;
zmax = 4;

rho = 3300;
alpha = 3e-5;

if exist('platform.s','file') == 0
  plat = 'native';
elseif exist('platform.s','file') == 2
  plat = 'ieee-le';
end
fid_t = fopen('time.0','r',plat);
fid1 = fopen('sxx.0','r',plat);
fid2 = fopen('syy.0','r',plat);
fid3 = fopen('temp.0','r',plat);
fid4 = fopen('aps.0','r',plat);
fid5 = fopen('phas.0','r',plat);

loadflac;
tsp = round(t(istp)/3600/24/365/1000);

for ii = 1:length(t)
    for jj = 1:nx
        %dZ = abs(diff([reshape(meshzz(1:end-1,jj,ii),ny,1);meshzz(1,jj,ii)+meshzz(end,1,1)]));
        dZ = -(diff([reshape(meshzz(1:end-1,jj,ii),ny,1);meshzz(1,jj,ii)+meshzz(end,1,1)]));
        drho = reshape(rho*(1-(alpha)*A4(:,jj,ii)),1,ny);
        tmass(jj) = drho*dZ*1e3;
    end
    zbox = meshzz(1,1,1)-meshzz(end,1,1);
    if tmass==0
        ziso(ii,:) = zeros(size(tmass));
    else
        ziso(ii,:) = zbox./(tmass./tmass(1))-zbox;
    end
end

%figure
clf
kt=istp;
a1=zeros(ny,nx); 

if istp > 3
    ma1 = (A1(:,:,kt-3) + A1(:,:,kt-2) + A1(:,:,kt-1) + A1(:,:,kt))./4;
    ma3 = (A3(:,:,kt-3) + A3(:,:,kt-2) + A3(:,:,kt-1) + A3(:,:,kt))./4;
    ma4 = (A4(:,:,kt-3) + A4(:,:,kt-2) + A4(:,:,kt-1) + A4(:,:,kt))./4;
else
    ma1=A1(:,:,kt);ma3=A3(:,:,kt);ma4=A4(:,:,kt);
end
ma5 = A5(:,:,kt);
xcnr=zeros(ny*nx,4); zcnr=xcnr;
xcnr(:)=Xcnr(:,:,kt);  zcnr(:)=Zcnr(:,:,kt);

vfac=4;qsize=1;
jv=[1:vfac:nx]; iv=[1:vfac:ny];

Vx=zeros([length(iv) length(jv)]); Vy=Vx; Xv=Vx; Yv=Vx;
Vx(:)=vx(iv,jv,kt); Vy(:)=vz(iv,jv,kt);
Xv(:)=meshxx(iv,jv,kt); Yv(:)=meshzz(iv,jv,kt);

subplot(311)
plot(meshxx(1,:,kt),meshzz(1,:,kt),'b','linewidth',2)
hold on
%plot(Xcnt(1,:,kt),ziso(kt,:),'k--')
ylabel('Topo (km)')
title(['Time = ',num2str(t(kt)/3600/24/365/1000),' kyrs'])
%axis([xmin xmax 0 1.5])
axis([xmin xmax -2 2])
[pos0,pos]=fixposition(1,-0.022,1,-0.05,.53,0,.6,0);
set(gca,'XTick',-45:15:45,'XTickLabel','')
%title('M = 0.9')

subplot(312)
patch(xcnr',zcnr',reshape(ma4',length(xcnr),1)','edgecolor','none')
hold on
quiver(Xv,Yv,Vx,Vy,qsize,'k');
%plot(meshxx(:,:,1),meshzz(:,:,1),'w.')
%contour(Xcnt(:,:,kt),Zcnt(:,:,kt),A4(:,:,kt),[600 1150],'k')
%plot([-100 100], [-6 -6],'color',[.8 .8 .8],'linewidth',2)
set(gca,'layer','top','box','on')
%caxis([0 0.5])
[pos0,pos]=fixposition(1,0,1,0,.5,0,1,0);
cb=colorbar;
cbxlb = get(cb,'Xlabel');
set(cbxlb,'string','APS')
cbpos=get(cb,'Position');
set(cb,'Position',[cbpos(1)+.12 cbpos(2) cbpos(3)*.5 cbpos(4)*.8])
ylabel('Depth (km)')
axis([-10 10 -3 0.6])
%axis([xmin xmax -25 5])
set(gca,'XTick',-45:15:45,'XTickLabel','')
set(gca,'Position',pos)
%for i=2:istp
%    hold on;
%    plot(.001*tracer(i,1),.001*tracer(i,2),'.') 
%end

subplot(313)
patch(xcnr',zcnr',reshape(ma3',length(xcnr),1)','edgecolor','none')
hold on
%contour(Xcnt(:,:,kt),Zcnt(:,:,kt),A4(:,:,kt),[600 1150],'k')
%plot([-100 100], [-6 -6],'color',[.8 .8 .8],'linewidth',2)
set(gca,'layer','top','box','on')
%caxis([0 1300])
[pos0,pos]=fixposition(1,0,1,0.05,.5,0,1,0);
cb=colorbar;
colormap jet;
cbxlb = get(cb,'Xlabel');
cbpos=get(cb,'Position');
set(cb,'Position',[cbpos(1)+.12 cbpos(2) cbpos(3)*.5 cbpos(4)*.8])
set(cbxlb,'string','Temp (^oC)')
xlabel('Across-Axis Distance (km)')
ylabel('Depth (km)')
axis([xmin xmax zmin zmax])
%axis([xmin xmax -25 5])
%set(gca,'XTick',-45:15:45)
set(gca,'Position',pos)
